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SYNCHRONOUS  SATELLITE  STATIONKEEPING  SIMULATION 


I.  INTRODUCTION 

The  motion  of  a  satellite  with  a  period  of  one  sidereal  day  is  perturbed 
by  the  gravitational  attraction  of  the  sun  and  moon,  radiation  pressure  and 
the  tangential  component  of  the  earth’s  gravitational  geopotential.  A  Syn¬ 
chronous  Satellite  Simulation  program  (SSS)  was  written  to  study  these  ef¬ 
fects  and  to  determine  the  feasibility  of  maintaining  a  satellite  near  a  speci¬ 
fied  synchronous  position  by  automatically  firing  rocket  motors  and  by  solar 
sailing.  The  different  parts  of  the  SSS  program  are  described  in  the  Appen- 
dice  s . 

II.  UNDERLYING  THEORY  TO  THE  COMPUTER  SIMULATION 

A.  Forces  Acting  On  the  Satellite 

The  influence  of  the  gravitational  attraction  of  the  sun  and  moon  contrib¬ 
ute  to  a  small  diurnal  oscillation  in  the  longitude  of  the  satellite  of  the  order 

of  0.  5  degrees  ^  at  the  most.  A  typical  daily  longitude  variation  has  been  cal- 

Z 

culated  by  Molitor  and  Kaplan  and  is  shown  in  Fig.  1.  However,  the  sun  and 

3 

moon’s  gravitational  attraction  does  cause  the  orbit  normal  to  precess  so 
that  the  inclination  changes  by  about  0.  85  degrees  per  year. 

Radiation  pressure  will  change  the  eccentricity  of  the  orbit  a  small 
amount  for  satellites  with  ordinary  area  to  mass  ratios  that  are  symmetric 
about  a  plane  that  contains  the  center  of  mass  of  the  satellite  and  is  perpen¬ 
dicular  to  the  equatorial  plane.  By  making  the  satellite  not  symmetric,  it  is 

possible  to  East- West  stationkeep  a  synchronous  satellite  by  radiation  pres- 

4 

sure  forces.  Methods  for  doing  this  are  presented  in  a  paper  by  Black  et  al  . 

B.  Orbital  Mechanics 

The  motion  of  a  satellite  is  determined  primarily  by  the  earth's  gravita¬ 
tional  potential  function 
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Fig.  1.  Typical  daily  longitude  variation  of  a  24-hour 
satellite  due  to  solar  and  lunar  perturbations. 
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For  the  purposes  of  studying  east-west  stationkeeping  the  effect  of  the 
gravitational  attraction  of  the  moon,  sun  and  planets  was  neglected.  The  SSS 
program  has  the  option  of  neglecting  terms  in  the  potential  function  greater 
than  order  two  or  three.  The  inaccuracies  in  measuring  and  producing  sur¬ 
faces  with  specified  reflective  properties  produce  uncertainties  in  the  forces 
acting  on  the  satellite  which  mask  the  above  effects. 

The  motion  of  a  24-hour  satellite  is  best  studied  by  using  an  earth  fixed 

5  7 

system  of  coordinates  rotating  at  the  earth's  angular  rate  U).  *  The  kinetic 

energy  of  the  satellite  expressed  in  this  rotating  system  is 

T  =  ~'m[r^4r^§4r^sin^0(A4U))^]  (2) 


Let  us  denote  the  forces  on  the  satellite  produced  by  rocket  thrusters  or  radi¬ 
ation  pre  s  sure  devices  as  F  ,  Fq,  F^.  Now  by  using  Lagrange’s  equations  and 
retaining  terms  to  the  second  order,  the  following  equations  of  motion  result. 


^  O  O  A 

r  -  r  8  -rsin  0  (  A  4  U))  =  -p/r^4  3fjiJ7Q  R*"  (3  cos  0-  l)/2r 

4  9  p  J  2  >  sin^  0  cos  2  (  A  -  A^^)/r^  +  F^/m  ( 3) 

3  n  J?nR2  3  J  RZ 

0  -  - - -  sin  2  0  -  - ^ -  sin  2  9  cos  2  (  A  -  A^ ^ ) 

2r  r 


4  —  ^  —  ■  sin  29  -  2r  9/r  4  F ,/mr  (4) 

. .  h  p  Jp  p  R  . 

A  =  - - -  sin2(A  -  A  ~ )  -  2  ctn  0  0  (  A  4  JO)  4  F  /mr  sin0  (5) 

D  L*L*  A 

r 

2r  (  A  4  :jj) 
r 

When  the  force  terms  are  set  to  zero  Fq,  F^,  and  F^  the  program  has  been 
found  to  simulate  the  motion  of  Syncom  II  to  a  high  degree  of  accuracy  as 
shown  in  Fig.  2.  The  discontinuity  in  the  curves  is  due  to  a  rocket  motor 
firing. 
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Fig.  2.  SYNCOM  II  drift 
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C.  Longitudinal  Perturbation  Force 

The  longitudinal  force  per  unit  mass  on  a  satellite  due  to  the  earth’s  geo¬ 
potential  is  given  by  taking  the  gradient  of  the  potential 
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For  an  equatorial  orbit  9  =  90  and 


-  8 

f.  /m  =  -  5.  568  x  10  sin  2  (X  -  X  ?)  newtons/kilogram 


22' 


(7) 


The  stable  equilibrium  points  are  located  at  the  positions  (X  -  X  ,^)  =  90°  and 
270°. 

If  a  satellite  were  placed  at  these  points,  the  satellite  wo\ild  remain  there 
or  oscillate  about  them  with  a  small  amplitude.  No  stationkeeping  apparatus 
would  be  required  if  the  satellite  were  placed  at  these  points  and  the  portion 
of  the  earth  "visible”  to  the  satellite  was  satisfactory.  The  stable  equilibrium 
points  are  located  off  the  western  coast  of  South  America  and  over  the  Indian 
Ocean.  However,  if  the  satellite  is  used  for  communications  between  the 
United  States  and  Europe  or  the  United  States  and  the  Far  East  the  satellite 
will  require  stationkeeping  apparatus.  The  visibilities  of  satellites  placed  for 
communication  between  the  United  States  and  Europe  and  between  the  United 
States  and  the  Far  East  are  shown  in  Fig.  3. 

III.  STATIONKEEPING  BY  SOLAR  SAILING 

The  instantaneous  radiation  pressure  forces,  on  the  sail  of  a  satellite, 
which  has  its  plane  in  the  north- south  direction  are  given  by  the  following  vec¬ 
tor  equations,  assuming  there  is  no  angular  dependence  of  the  solar  absorp¬ 
tivity.  These  equations  hold  for  any  orbit. 

Fx  =  T,h,2-a),v  3)2 1 -sV  ?i 

-  sOt  |  SN  •  U  |  (U-  S)  sin  [cos'1  (-S^  •  ?)])  (8) 
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of  satellites  placed  at  stable  equilibrium  points. 
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where  h,  s,  g,  f  are  +1  with  the  following  exceptions: 
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A.  Sail  Type  I 

The  simplest  type  of  sail  for  a  satellite  that  has  one  side  oriented  to  the 
earth  and  one  side  oriented  toward  the  north  direction,  is  a  sail  whose  plane 
is  coincident  with  the  radius  to  the  satellite  from  the  center  of  the  earth  and 
the  north- south  direction.  (See  Fig.  4)  This  sail  has  a  low  solar  absorptiv¬ 
ity  on  one  side  and  a  high  solar  absorptivity  on  the  other.  The  thermal  emis- 
sivity  of  both  sides  is  approximately  the  same.  Under  these  conditions  the 
force  resultant  from  thermal  radiation  from  both  sides  of  the  sail  cancels  for 
a  thin  sail.  The  net  force  along  the  orbit  is  then  due  only  to  solar  radiation 
pressure.  The  average  annual  force  in  the  longitudinal  direction  is  given  by 


Sa 

c 


<  cos 


> 


(11) 


where  <  cos  i  >  is  average  annual  effect  of  the  sunTs  movement  above  and 
below  the  orbital  plane  of  the  satellite.  If  one  assumes  that  the  earth’s  or¬ 
bit  about  the  sun  is  circular  and  an  equatorial  orbit  for  the  satellite  <cos  i>  « 
0.  96  and  3^  is  the  solar  ab sorptivity  of  the  sail  facing  in  the  direction  of  orbital 
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Fig.  4.  Solar  sail  Type  I. 
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motion  and  is  solar  absorptivity  of  the  opposite  face.  This  type  of  sail  does 
not  lend  itself  to  being  adapted  to  the  stationkeeping  methods  to  be  discussed 
in  Section  IV. 

B.  Sail  Type  II 

The  second  type  of  sail  is  similar  to  the  first  except  both  sides  of  the  sail 
have  the  same  highly  reflective  surface  and  can  be  rotated  along  the  earth  sat¬ 
ellite  line  so  that  the  plane  of  the  sail  goes  from  the  north- south  to  east-west 
direction.  The  sail  is  left  in  the  north- south  portion  for  half  the  orbit  starting 
(or  ending)  with  the  sun-earth  line  and  the  east-west  position  for  the  other  half 
of  the  orbit.  (See  Fig.  5)  The  average  annual  force  in  the  longitude  direction 
is 

F  ,  =  I"  ±  —  — )  <  cos  i  >  -  a  <  sin  2i  >  1  (12) 

aX  L  4  c  Zttc  J  v 

The  positive  sign  is  used  when  the  north- south  position  is  taken  for  half  the  or¬ 
bit  starting  with  the  earth  sun  line  and  the  negative  sign  is  used  when  the  north- 
south  position  ends  with  the  earth  sun  line.  This  type  of  sail  is  more  compli¬ 
cated  than  the  first  type  because  it  requires  a  rotating  mechanism  and  sensor 
logic  to  determine  the  location  of  the  earth  sun  line.  It  does  require  less  sail 
area  than  the  first  type.  It  is  also  adaptable  to  the  previously  mentioned  meth¬ 
ods  of  thrusting  by  rotating  the  sail  into  or  out  of  position  with  respect  to  the 
earth-  sun  line . 

C.  Sail  Type  III 

The  third  type  of  sail  has  both  sides  highly  reflective.  It  rotates  about  the 
north- south  direction  continuously  with  a  period  of  12  hours.  (See  Fig.  6)  The 
average  annual  force  in  longitudinal  direction  is 

2 

F  =  —  -r  +  Z  (1  -  a)  sin  2^p  <cos  i  >  (13) 

where  cp  is  the  angle  between  the  plane  of  the  sail  and  the  projection  of  the  earth 
sun  line  on  the  equatorial  plane  when  the  satellite  is  crossing  the  projection  of 
the  earth  sun  line  on  the  equatorial  plane.  This  type  of  sail  requires  the  least 
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Fig.  5.  Solar  sail  Type  II. 


Fig.  6.  Solar  sail  Type  III. 
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amount  of  sail  area  for  stationkeeping.  If  one  sets  CP  =  45°  and  assumes  a  solar 
absorptivity  of  0.  9  for  one  side  of  the  Type  I  sail  and  0.  05  for  the  other  types 
of  sails,  the  relative  area  needed  to  produce  the  same  average  annual  force  in 
the  longitudinal  direction  for  Type  I,  Type  II  and  Type  III  sails  is  in  the  ratio 
of  3.  84  :  1.  70  :  1  respectively.  A  typical  satellite  with  a  mass  of  100  Kg  would 
require  a  Type  III  sail  with  an  area  of  1.47  square  meters  to  stationkeep  at  the 
position  of  highest  longitudinal  force  due  to  the  earth’s  geopotential.  By  varying 
the  value  of  cp  it  is  possible  to  vary  the  longitudinal  force  along  or  against  the 
orbital  motion.  Therefore,  this  type  of  sail  is  adaptable  to  any  method  of  auto¬ 
matic  stationkeeping. 

D.  Results  of  Computer  Simulations 

The  computer  simulation  of  stationkeeping  a  synchronous  satellite  by  solar 
sailing  revealed  little  difference  between  the  three  types  of  sail  considered  for 
stationkeeping.  A  typical  stationkeeping  simulation  is  shown  in  Fig.  7.  In  this 
simulation  sail  Type  III  was  used  with  =  -45°  and  a  =  .  1.  The  area  of  the  sail 
was  1.35  square  meters  for  a  100  Kg  satellite.  The  satellite  was  stationkept 
near  a  position  at  which  a  maximum  longitudinal  force  occurs.  (296°  E  longitude) 

IV.  AUTOMATIC  STATIONKEEPING 

A.  Stationkeeping  by  Constant  Force 

Different  methods  of  applying  a  tangential  force  to  counteract  the  longitudi¬ 
nal  perturbational  force  were  simulated  by  the  program.  The  simplest  of  these 
is  to  fire  a  rocket  motor  with  the  same  value  of  longitudinal  force  that  is  acting 
upon  the  satellite.  However,  in  practice  the  force  cannot  be  exactly  matched 
for  each  longitudinal  position  of  the  satellite.  The  method  is  limited  if  there  is 
no  control  over  the  force  level,  once  the  satellite  is  placed  in  a  longitudinal  po¬ 
sition  within  the  longitude  regions  45°  <  X  -  <  135°  and  225°  <  X  -  <315°. 

If  the  force  level  is  not  exactly  matched  for  a  particular  longitudinal  position  the 
satellite  will  undergo  an  oscillation  about  the  longitude  for  which  it  is  exactly 
matched.  A  simulation  of  an  unmatched  force  level  was  done  by  the  SSS  program 
and  .is  shown  in  Fig.  8.  The  satellite  was  initially  at  291°  longitude  and  the  force 
level  was  equal  to  the  longitudinal  force  at  283.  5°  longitude. 
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Fig.  7.  Variations  in  semi-major  axis 
and  longitude  versus  time. 
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Fig.  8.  Simulation  of  an  unmatched  force  level. 
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B.  Automatic  Stationkeeping  -  Unidirectional  Firing 

g 

Automatic  stationkeeping  was  proposed  by  W.  E.  Morrow  in  a  memo  dated 
1  July  1965.  A  block  diagram  of  the  logic  required  by  the  satellite  for  intern¬ 
ally  computing  its  position  and  rocket  firing  times  is  shown  in  Fig.  9.  This 
description  does  not  describe  the  LES-6  stationkeeping  logic. 

It  is  the  purpose  of  the  coincidence  logic  to  provide  a  pulse  derived  from 
fan  beam  sensors  at  a  predetermined  satellite  position  with  respect  to  the  earth 
and  sun.  The  satellite  may  be  assumed  spinning  so  that  the  output  from  the 
sensors  are  in  the  form  of  pulses  or  for  non- spinning  satellite  the  output  from 
the  sensors  may  be  chopped  to  produce  the  pulses.  This  coincidence  pulse  per¬ 
mits  comparison  between  a  running  clock  and  a  fixed  clock  in  the  decision  logic. 
Almost  any  position  in  the  orbit  can  be  used  for  establishing  the  satellite’s  po¬ 
sition  where  the  optical  sensors  can  view  an  illuminated  portion  of  the  earth. 
For  definiteness  let  us  pick  the  6  a.m.  position  shown  in  Fig.  10. 

The  decision  logic  determines  whether  the  satellite  is  early  or  late  in  ar¬ 
riving  at  the  M6  a.m.  position”  each  day.  The  device  compares  the  coincidence 
signal  with  the  output  signal  from  the  clock  logic. 

If  the  signal  from  the  coincidence  logic  occurs  in  position  A  of  Fig.  1  1,  the 
satellite  is  early  for  the  three  days  shown  and  no  signal  is  sent  to  the  thrust 
logic.  If  a  signal  from  the  coincidence  logic  falls  in  position  B  then  a  signal  is 
sent  to  the  thrust  logic. 

The  thrust  logic’s  purpose  is  to  fire  the  rocket  motor  in  the  direction  of 
the  orbit  tangent  for  a  specified  length  of  time  near  the  6  a.m.  position.  This 
is  repeated  daily  until  the  satellite  is  "on  station”.  The  thrusting  is  then 
stopped  until  another  signal  is  received. 

The  clock  must  have  an  accuracy  of  5  in  10^/day  for  5  years  and  put  out 
a  signal  (alarm)  at  10  hr  1Z  min  GMT  each  day.  This  signal  is  sent  to  the  de¬ 
cision  logic  for  further  processing. 

The  purpose  of  the  clock  logic  is  to  apply  a  daily  correction  to  the  clock 
to  allow  for  the  fact  that  the  sun  is  ahead  of  or  behind  the  clock  according  to 
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Fig.  9.  Satellite  logic  for  computing  position  and  rocket  firing  times 
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Fig.  10.  6a.m.  position  of  satellite  in  orbit. 
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Fig.  11.  Signal  from  clock  logic. 


16 


the  time  of  year.  The  correction  that  must  be  added  to  the  time  is  shown  in 
Fig.  12.  The  curve  shown  in  Fig.  12  has  to  be  approximated  by  a  function  that 
is  easy  to  program  into  the  satellite  logic.  The  errors  associated  with  the 
degree  of  exactness  to  which  the  clock  logic  will  correct  for  the  difference  be¬ 
tween  apparent  solar  and  solar  mean  time  were  accounted  for  from  tables  fur¬ 
nished  by  A.  Braga-Illa  for  LES-6. 

As  an  example  of  how  the  correction  is  applied,  let  us  consider  the  cor¬ 
rection  at  1  5  January  1967.  To  do  an  exact  calculation  of  the  time  to  be  added 
we  would  have  to  go  through  an  iterative  process.  However,  for  our  purposes 
the  correction  may  be  taken  as  the  correction  that  occurs  on  January  1  5.  0  which 
is  9.  0  min.  Therefore,  the  clock  logic  is  to  put  out  a  signal  to  the  decision 
logic  which  starts  at  10  hr  21  min  GMT  and  lasts  for  one  hour.  When  the  sun 
is  ahead  of  the  clock  the  correction  has  to  be  subtracted  from  24  hours  and  the 
result  added  to  the  time  of  the  previous  days  output. 

The  method  is  a  little  cumbersome  in  logic  design.  A  better  way  to  do  it 
would  be  to  shift  0  time  base  line  so  that  all  corrections  become  additions. 

It  is  possible  to  do  without  the  clock  logic  by  using  an  optical  sensor  that 
senses  the  declination  of  the  sun  and  relative  position  of  the  earth  satellite  line 
at  the  same  time  to  correct  for  the  actual  longitudinal  motion  of  the  sun  rela¬ 
tive  to  the  "mean"  sun  in  the  celestial  sphere.  The  sensor  can  be  visualized 

as  having  a  cluster  of  pencil  beams  arranged  in  the  familiar  analemma  found 

9 

on  some  globes.  A  sensor  of  this  type  is  being  constructed  by  C.  Burrowes 
for  use  in  LES-6. 

C.  Effects  of  Sensor  and  Clock  Errors 

There  are  errors  associated  with  the  measurement  of  longitude  due  to  the 
beamwidth  of  the  sensors  and  the  degree  of  exactness  to  which  the  clock  logic 
will  correct  for  the  difference  between  apparent  solar  and  solar  mean  time. 
There  are  also  errors  induced  because  the  LES-6  satellite  spin  axis  is  not 
aligned  perfectly  in  the  north- south  direction.  The  source  of  error  has  been 
investigated  by  B.  J.  Moriarty.  ^  The  shape  of  the  probability  density  func¬ 
tion  of  the  sum  of  the  errors  is  not  known.  Therefore,  a  uniform  probability 
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Fig.  12.  Ephemeris  correction  vs.  time  of  the  year. 
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function  with  a  width  of  one  degree  of  longitude  was  assumed  for  use  in  the  SSS 
program.  The  errors  associated  with  the  analemma  sensor  for  LES-6  were 
separately  accounted  for  from  tables  furnished  by  C.  Burrowes. 

The  influence  that  the  errors  have  on  the  motion  of  the  satellite  is  shown 
in  Figs.  13  and  14.  Fig.  13  is  the  result  of  an  SSS  computer  simulation  of  an 
automatic  stationkeeping  system  using  unidirectional  firing  with  no  errors  pres¬ 
ent  in  determining  longitude.  Fig.  14  is  the  result  of  a  simulation  under  the 
same  conditions  as  the  simulation  of  Fig.  13  except  the  clock  logic  errors  and 
sensor  errors  are  introduced.  The  errors  do  not  significantly  affect  the  range 
of  longitudes  over  which  the  satellite  may  be  stationkept.  The  errors  in  longi¬ 
tude  sensing  were  less  than  0.  7  degrees  due  to  clock  errors  and  less  than  0.  5 
degrees  due  to  random  sensor  errors. 

D.  Automatic  Stationkeeping  Unidirectional  Firing  With  Force  Cut-Down 

The  best  method  of  automatic  stationkeeping  so  far  proposed  is  that  pro¬ 
posed  by  Roger  Brockett.^  This  method  uses  the  least  amount  of  fuel  and  sta- 
tionkeeps  within  the  smallest  range  of  longitude.  As  an  example  of  how  this 
method  operates,  let  us  assume  the  satellite  is  initially  located  at  a  longitude 
such  that  the  earth’s  gravitational  component  along  the  orbit  will  cause  the  lon¬ 
gitude  of  the  satellite  to  decrease  in  time.  The  longitude  of  the  satellite  is  al¬ 
lowed  to  decrease  until  it  goes  past  a  selected  longitude,  X  .  At  this  point  the 
rocket  motors  are  fired  for  a  specified  length  of  time  such  that  the  impulse  given 
per  day  is  approximately  five  times  the  "impulse”  given  to  the  satellite  by  the 
earth's  geopotential.  This  is  the  same  as  the  method  previously  described. 
However,  the  longitude  is  sampled  once  per  day  and  its  value  is  stored  in  a 
memory  bank  in  the  satellite.  If  the  longitude  measured  on  a  given  day  is 
greater  than  that  measured  the  day  before,  the  rocket  impulse  given  per  day 
is  reduced  until  the  satellite  again  passes  the  selected  longitude,  X  ,  at  which 
time  the  thrusting  is  stopped.  The  amount  that  the  impulse  may  be  reduced  is 
a  variable  that  may  be  put  in  the  input  to  the  SSS  program.  A  typical  example 
of  this  stationkeeping  method  is  shown  in  Fig.  15. 
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Fig.  13.  Results  of  automatic  stationkeeping  with  unidirectional 
firing  with  no  errors  in  determining  longitude. 
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Fig.  14.  Results  of  automatic  stationkeeping  with  unidirectional 
firing  with  clock  logic  error  and  sensor  errors  introduced. 
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Fig.  15.  Results  of  automatic  stationkeeping  with  unidirectional 
firing  with  force  cut-down. 
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E.  Automatic  Stationkeeping  Bidirectional  Thrusting 

It  is  possible  to  stationkeep  with  a  system  which  thrusts  in  one  direction 
when  the  satellite  exceeds  a  given  longitude  and  in  the  other  direction  when 
the  satellite  is  less  than  another  given  longitude.  In  other  words,  whenever 
the  satellite  has  drifted  out  of  a  given  band  of  longitude,  thrusting  is  done  to 
return  the  satellite  to  the  band.  This  method  of  stationkeeping  is  not  recom¬ 
mended  for  use  with  thrust  systems  that  consume  fuel  at  a  low  specific  im¬ 
pulse  because  a  lot  more  fuel  is  needed  than  with  the  methods  previously  men¬ 
tioned.  However,  it  is  adaptable  to  solar  sailing  techniques  where  no  fuel  is 
used.  A  simulation  of  the  bidirectional  technique  is  shown  in  Fig.  16.  This 
simulation  was  run  under  the  same  conditions  as  the  unidirectional  firing  si¬ 
mulation  of  Fig.  13.  The  fuel  consumption  for  stationkeeping  a  132  Kg  satel¬ 
lite  for  1000  days  was  3.  85  Kg  for  a  rocket  system  which  has  an  I  =70 
for  bidirectional  thrusting.  Under  the  same  conditions  only  0.86  Kg  are  needed 
for  unidirectional  thrusting. 

Another  method  of  bidirectional  thrusting  has  been  proposed  by  A.  Braga- 
Illa.  The  logic  of  this  method  is  explained  in  Appendix  E  describing  subrou¬ 
tine  Newfor  of  the  SSS  program.  It  is  well  suited  for  solar  sailing  techniques. 

A  simulation  of  this  method  is  shown  in  Fig.  17. 
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Fig.  1  6.  Results  of  automatic  stationkeeping  using  bidirectional  firing. 
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Fig.  17.  Results  of  automatic  stationkeeping  using  A.  Braga-Illa's 
force  scheme  . 
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APPENDIX  A 


MAIN  PROGRAM  OF  (SSS) 


The  main  routine  of  the  SSS  program  issues  calls  to  the  subroutines, 
solves  the  equations  of  motion,  and  produces  the  output  of  the  program 
(Fig.  A- 1  and  Appendix  J).  The  initial  section  of  comment  cards  gives  a 
brief  description  of  the  program,  its  control  variables  and  howto  run  it. 

The  first  executable  statements  are  definitions  of  constants.  The  Main  pro¬ 
gram  calls  GETINP  to  obtain  the  input  data  and  control  variables.  After  a 
sequence  of  housekeeping  details  --  moving  variables  to  and  from  commons, 
changing  variable  types,  and  changing  variables  in  degrees  to  radians  --  the 
values  necessary  for  the  first  entry  into  NEWFOR,  FLAM,  and  CIND  are 
initialized.  These  variables  are  initialized  after  the  call  to  GETINP  so  that 
they  are  reinitialized  at  the  beginning  of  subsequent  sets  of  input  data.  The 
control  variable  II  determines  if  the  branch  using  second  harmonics  in  the 
equations  of  motion  or  the  branch  using  third  harmonics  is  taken.  In  either 
case  calls  are  issued  to  SSL  and  FLAM.  SSL  initializes  CIND  for  solution 
of  differential  equations;  FLAM  returns  the  force  in  the  A  direction.  Certain 
logical  paths  determined  by  the  input  data  return  other  variables  to  Main 
through  COMMON  (see  Appendix  D).  After  the  call  to  FLAM,  D2R,  D2XLAM 
and  D2THET  are  calculated.  Both  branches  meet  at  statement  23,  where 
three  checks  are  made.  First,  is  this  time  through  the  end  of  a  day?  Sec¬ 
ond,  is  this  time  through  the  end  of  ten  days?  Third,  is  this  time  through 
the  end  of  amount  of  time  asked  for?  If  answer  to  all  three  questions  are  no, 
then  it  is  end  of  one  of  the  86,  400  times  per  day  that  the  variables  D2R, 
D2XLAM,  and  D2THET  are  calculated.  At  statement  6  the  calls  to  FINDV, 
and  DPNV  provide  T,  R,  THETA,  XLAM,  DR,  DTHETA,  and  DXLAM  for 
the  following  loop.  Again  a  branch  is  taken  based  on  II.  If  it  is  the  end  of 
a  day,  the  following  are  printed  out. 


R 


DR 


TDEG 


TDAYS 


days  since  beginning  of  this  set  of  input  data 
radius  of  orbit 

change  in  radius  in  last  loop  of  the  day 
THETA  in  degrees 


26 


DTDEG 

XLDEG 

DXLDEG 

XELONG 


change  in  THETA  in  last  loop  of  the  day 
XLAM  in  degrees 

change  in  XLAM  in  last  loop  of  the  day 
east  longitude  in  degrees 


After  this  printout  a  branch  is  made  to  question  three.  If  the  answer  to  ques¬ 
tion  two  is  yes  also,  there  precedes  this  printout  a  line  of  print  containing  the 
following  averaged  over  previous  10  days. 


A 

radius 

E 

eccentricity  of  the  orbit 

XI 

inclination 

XIM1D 

minimum  inclination 

XLMAD 

maximum  inclination 

DIFFER 

XLMAD  -  XLMID 

Every  ten  days  there  is  also  punched  a  card  containing  TDAY,  A,  and  XLAM, 
to  be  used  by  a  plotting  program.  After  this  there  is  a  branch  to  question  3. 

If  the  answer  to  question  3  is  yes,  and  the  amount  of  output  for  one  set  of  data 
is  finished,  the  program  branches  to  7,  CALL  GETINP.  The  last  set  of  data 
contains  NUM  =  -1,  which  causes  program  to  terminate  normally.  Each  set 
of  data,  except  the  last,  takes  11/2-2  l/2  hours  of  computer  time. 


27 


K/l  MN 


PK.OSK.AM 


f="l_OW 


28 


APPENDIX  B 


The  subroutine  GETINP  establishes  the  logical  variables,  all  the  input, 
and  the  modes  of  output.  There  are  eight  classes  of  integers  that  are  Tised  in 
conjunction  with  the  logical  variables:  1 1 ,  1 2,  13,  14,  15,  16,  17,  18.  All  the 
parameters  are  initialized  in  GETINP  and  then  those  that  are  to  be  changed 
are  read  in  through  a  namelist  setups.  Using  the  logical  parameters  and  the 
numerical  parameters ,  the  subroutine  SETUP  and  its  entry  point  RITOUT 
prepare  the  variable  format  for  PRINT  and/or  PUNCH  so  that  at  the  begin¬ 
ning  of  each  output  one  may  quickly  review  the  particular  set  of  parameters 
used.  Control  is  then  returned  to  the  main  program  which  then  prints  the 
titles  for  the  output,  prints  and/or  punches  the  orbital  information  every  ten 
days . 

11  determines  whether  the  calculations  in  the  main  program  use  second 
order  harmonics  or  third  order  harmonics.  11  is  set  to  1  which  selects  the 
second  order  harmonics.  If  third  is  true,  12  then  equals  2  and  the  third  order 
harmonics  are  used.  12  controls  the  punched  output.  12  is  initially  1  which 
produces  punched  output.  When  no  punched  output  is  desired,  12  should  read 
in  as  2.  The  punched  output  is  used  by  a  separate  program  Z-PLOT  which 
produces  one  frame  of  alphanumeric  data  and  a  second  frame  with  2  graphs 
longitude  vs.  days  and  semi-major  axis  vs.  days.  13  controls  the  solar  time 
correction  and  primitive  uses  of  a  solar  sail.  13  is  initialized  as  1,  no  cor¬ 
rection.  If  13  equals  2  or  linear  is  true,  there  is  a  linear  solar  time  correc¬ 
tion.  If  13  equals  3  or  LINRAN  is  true,  there  is  a  linear  solar  time  correc¬ 
tion  with  random  error.  If  13  equals  4  or  perfect  is  true,  the  solar  time  is 
corrected  perfectly.  If  13  equals  5  or  STEP  is  true,  the  solar  time  is  corrected 
by  a  step-function.  When  13  equals  6  or  STPRAN  is  true,  the  step  function 
correction  is  modified  by  calling  the  random  error  function.  If  13  equals  7  or 
EPH  is  true,  the  solar  time  correction  is  done  by  an  ephemeris  sun  sensor. 
When  13  equals  8  or  EPHRAN  is  true,  the  ephemeris  sun  sensor  correction  is 
adjusted  using  the  random  error  function.  The  remaining  13  values  pertain  to 
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the  solar  sail.  When  13  equals  9  or  SAIL1  is  true,  a  fixed  sail  is  used  with  a 
constant  force  along  lambda.  If  13  equals  10  or  SAIL3  is  true,  a  fixed  sail 
scheme  3  is  used.  14  determines  which  XJ22  and  XLAM22  are  to  be  used. 

The  program  is  set  up  to  use  14  equal  to  2  which  sets  XJ22  at  -  1.700E-6  and 
XLAM  at  -19  degrees.  There  are  no  logical  variables  that  go  with  14.  15  es¬ 

tablishes  the  force  cutdown.  Initially  15  equals  1,  no  force  cutdown.  When  15 
equals  2,  there  is  a  force  cutdown  equal  to  the  input  number  CUT.  To  obtain 

15  equals  2,  it  can  be  read  in  directly  or  CUTDOWN  can  be  read  in  as  true. 

16  and  18  are  the  more  sophisticated  uses  of  the  solar  sail.  16  equal  to  1  or 

SAIL2  equal  to  true  fixes  the  sail  and  allows  one  to  also  make  a  sun  correc¬ 
tion.  16  is  initialized  as  0.  16  equal  to  2  or  TACSAL  is  true  the  double  pre¬ 

cision  function  FLAM  calls  TACK  and  the  tacked  sail  scheme  is  used.  18 
flips  the  sail  according  to  three  schemes:  18  equal  1  or  FLIPW  equal  to  true 
fixes  the  sail  through  half  the  orbit  and  flips  it  off  with  the  orbit;  18  equals  2 
or  FLIPAG  true  fixes  the  sail  through  half  the  orbit  and  flips  if  off  against  the 
orbit;  18  equals  3  or  SAIL4  true  flips  the  sail  due  to  logic  within  the  satellite. 

17  determines  if  the  orbit  elements  are  averaged  at  the  end  of  every  ten  days 
before  they  are  punched.  Initially  17  equals  2  and  AVERAG  is  false,  which 
loops  the  main  program  around  the  averaging  sequence.  When  17  equals  1  or 
AVERAG  is  true,  the  main  program  averages  every  ten  days  and  punches  the 
averaged  figures. 

The  basic  units  for  the  parameters  are  kilometers,  seconds,  degrees. 
They  are  initialized  as  follows: 


A 

A1 

A2 

ALONGR  = 
CONST  = 


42165.  5D0 
0.  9D0 
0.  1D0 
0.  0 

-4.  45D-8 


radius  of  earth  in  kilometers 
solar  absorptivity  for  solar  sail 
solar  abosrptivity  for  solar  sail 
force  along  the  radius  to  the  satellite 
kilonewtons/kilogram 

=  1.  6  *  (factor  =  5.  6)/365.  2  5/(86.  4  *  6) 

sign  is  negative  in  second  and  fourth 
coordinates 
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CUT 

— 

0.  0 

DATE 

0.  0 

E 

= 

o.  0 

FORCCT 

= 

0.  0 

force  along  theta 

FORLAM 

1.  0 

force  along  lambda 

IRUN 

— 

0 

a  counter  so  one  can  keep  track  of  the 
various  computer  runs 

ITHDAY 

= 

10 

frequency  in  days  of  print  and/or  punch 
execution 

NUM 

6 

number  of  impulses  given  (used  in  FLAM) 

OFFSET 

-90.  0D0 

angle  by  which  sail  can  be  offset  in  the 
TACK  scheme 

OL 

= 

0.  0D0 

argument  of  perigee  for  first  orbit,  in 
degrees 

TEND 

— 

1800.  0D0 

number  of  solar  days  program  is  to  run 
for 

XDISP 

70.  0D0 

specific  impulse  of  thruster  in  seconds 

XI 

= 

0.  0D0 

inclination  of  the  orbit  in  degrees 

XLAMD 

=: 

300.  0D0 

degrees  of  longitude  measure  east  of 
Greenwich 

XLIMD 

297.  0D0 

minimum  angle  for  thrusting,  simple 
scheme 

XLIMXD 

— 

350.  0D0 

maximum  angle  for  thrusting,  simple 
scheme 

XMASS 

= 

132. 0D0 

mass  of  satellite  in  kilograms 

In  addition  the 

re 

is  a  logical  parameter,  NEWFR,  that  has  no  numerical  equi- 

valent.  NEWFR  is  initially  FALSE.  When  NEWER  is  true,  the  SUBROUTINE 
NEWFOR  is  called  by  FLAM  and  the  magnitude  of  the  impulse  if  modified  ac¬ 
cordingly. 

To  change  any  of  the  parameters  or  logical  variables  the  NAMELIST/ 
SETUPS/  is  used  &  SETUPS  is  on  the  first  card  beginning  of  column  2.  The 
variables  desired  changes  are  then  punched  separated  by  commas.  To  set  a 
logical  variable  to  TRUE  just  use  T  for  example:  NEWFR  =  T  would  be 
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sufficient  to  set  NEWFR  true.  Do  not  punch  in  column  73  or  beyond.  When 
using  a  second  card  begin  in  column  3.  Do  not  follow  the  last  variable  with  a 
comma;  instead  leave  one  or  more  spaces  and  punch  &  END.  The  set  of  cards 
with  the  SETUPS  changes  is  then  placed  behind  the  //FT0  5F001  control  card. 

The  parameters  are  printed  and/or  punched  according  to  a  variable  FOR¬ 
MAT  routine,  RITOUT.  RITOUT  is  an  ENTRY  point  in  SETOUT.  The  first 
call  to  SETOUT  establishes  whether  or  not  to  punch  when  RITOUT  is  called. 
RITOUT  is  called  for  each  I  variable.  To  add  an  I  variable  to  the  program, 

IY,  where  Y  is  an  integer  greater  than  8,  eight  additions  to  GETINP  are  nec¬ 
essary.  First  a  comment  card  to  describe  the  IY  variable  and  the  X  uses  to 
which  IY  can  be  put.  Secondly,  it  must  be  initialized.  The  remaining  additions 
are  to  setup  the  call  to  RITOUT.  There  are  X  +  1  dimension  statements  as 
follows:  Dimension  VFIY(16,6)  and  X  statement  like  dimension  XVFIYX(l6). 

For  each  X  there  is  also  an  EQUIVALENCE  statement  such  as:  EQUIVALENCE 
(XVFIXY(l),  VFIX(1,  Y).  For  each  value  IX  can  assume  there  must  be  a  data 
statement  in  Hollorith  form  containing  the  alphanumeric  data  you  desired 
printed  and/or  punched  to  describe  the  alternative  used.  You  need  to  add  IX 
to  the  101  format  statement,  and  one  call  to  RITOUT,  for  example:  Call  RIT¬ 
OUT  ( VFIY ( 1 ,  X)) . 

GETINP  thus  establishes  the  particular  parameters  and  logical  choices 
for  one  computer  run,  which  takes  approximately  two  hours  on  the  IBM/360 
with  TEND  =  1800. 
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APPENDIX  C 


INTEGRATION  FUNCTION  AND  SUBROUTINE 


The  double  precision  function  CIND(K,  A,  B)  and  its  related  functions 
FINDV,  DPNV,  and  SUBROUTINE  SSL  find  the  first  and  second  derivatives 
and  calculate  new  values  for  R,  THETA,  XLAM,  DR,  DTHETA,  DXLAM, 
and  T  for  each  cycle  through  the  main  program.  The  call  to  SSL  sets  up  the 
constants  for  SSL.  FIND  calls  CIND  with  K  =  0  which  routes  the  control  to 
the  proper  place  to  find  a  new  T.  For  the  other  variables  DPNV  is  used. 

In  the  computer  GO  TO  statements  only  the  first  two  solutions  are  used  by 
this  program.  The  sequence  of  four  subprograms  was  written  in  another 
language  outside  Lincoln,  so  there  are  more  choices  than  this  program 
uses. 

The  integration  formula  used  is  the  Adams  four-point  formula: 

yn+i  =  yn  +  ik  <55^  -  59Ci +  37Cz  -  9^-3>  <c-‘> 

where  h  =  x  ^  -  xn  is  the  constant  incremental  change  of  the  independent  vari¬ 
able  x,  y  being  a  function  of  x. 

Starting  formulas  are  needed  to  obtain  the  first  three  points.  The  follow¬ 
ing  set  of  formulas  is  used  to  obtain  y ^  y^,  and  y^: 


vi1' 

=  yo  +  hV 

(C-2) 

yi2) 

=  +!  <y(i1)' +  yo,} 

(C-  3) 

vi3) 

=  y0  +z  (2y(i2)'  +y(i1)' +  3y0,) 

(C-4) 

=  y(i3)  + 1  (^'i31,  -  y0'> 

(C-  5) 
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(C-  6) 


y 


(i) 

3 


y 


(2) 

3 


yi3,+4(5y(.,.+8y(3),.yo,) 

42)  +  i!  <23yf'-  S3)' +  %'> 

y(2Z)  +  Ta  {9Y(31]'  +  19^z2)'  -  5^i3)' +  V} 


(C-  7) 
(C-  8) 


where  the  superscripts  indicate  the  order  of  the  approximation.  The  past  val¬ 
ues  of  the  derivatives  are  stored  in  reserved  blocks.  For  a  more  complete 
discussion  of  the  starting  equations,  the  reader  is  invited  to  contact  the  author 
In  particular,  a  substitution  of  equations  (1)  and  (2)  and  similar  expressions 
for  their  derivatives  into  equation  (3)  will  show  that  equation  (3)  reduces  to  the 
first  four  terms  of  the  Taylor  series. 

The  truncation  error  of  formula  I  is  tLttt  h  (s),  where  x  <  s  x  ,  ,  . 

1 2  0  y  n  n  + 1 

The  accumulation  of  the  increments  (Ay)  =  4^-  (S5y  1  -  39y!  ,  +  37y!  . 

n  24  7  n  n-  1  n-2 

-  9y!  is  performed  in  double  precision  to  substantially  reduce?  the  round-off 
error. 

Complete  unde rflow- ove rflow  testing  is  made.  An  underflow  is  corrected 
and  a  zero  stored. 
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DOUBLE  PRECISION  FUNCTION  CIND  (K.,A,B) 
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APPENDIX  D 


DOUBLE  PRECISION  FUNCTION  FLAM  (XLAM) 


The  double  precision  function  FLAM(XLAM)  computes  a  force,  FLAM, 
simulating  the  force  of  the  thruster  or  solar  sail  depending  on  the  logic  with 
which  the  program  is  run.  FLAM  is  returned  to  the  main  program.  The 
force  is  applied  NUM  times  an  orbit,  controlled  by  the  gate  K.  When  K  is 
positive  the  gate  is  closed  returning  both  force  and  FLAM  as  0.  When  K  is 
negative  control  passes  to  label  8.  The  various  routes  for  the  logic  depend 
on  the  force  scheme  used  for  the  particular  run  in  question.  The  variable 
15  controls  the  cutdown;  when  15  is  1  there  is  no  cutdown.  When  15  is  2,  the 
force  is  cut  by  the  factor  CUT,  as  long  as  LLL  is  also  positive.  LLL  is  the 
control  variable  whose  sign  is  determined  by  the  satellites  position  within  a 
specified  band,  delineated  by  XLIM  and  XLIMX.  LLL  is  negative  when 
XLAM  is  inside  the  band  and  positive  when  XLAM  is  outside  the  band.  13, 
for  values  1-8,  controls  the  type  of  clock  error  returned  by  CORR(TEL). 
From  this  an  XTEMP  is  calculated.  XLAM  is  the  actual  position  of  the  sat¬ 
ellite;  XTEMP  is  where  the  satellite’s  logic  thinks  it  is.  The  old  force  cal¬ 
culations  use  LLL,  whereas  the  new  force  calculations  use  the  subroutine 
NEWFOR.  As  long  as  the  counter  M  is  less  than  NUM,  control  returns  to 
the  main  program  with  FLAM  =  FORCE.  Once  an  orbit  the  fuel  consumption, 
the  total  number  of  impulses,  JFK,  CORR  the  corrected  longitude  in  degrees 
are  printed  out.  The  fuel  consumption  is  kept  as  a  running  sum.  After  the 
fuel  report  is  printed,  on  the  last  loop  during  which  there  is  thrusting,  K  is 
set  positive  and  the  other  counters  are  initialized  for  the  next  orbit. 

For  the  various  schemes  of  solar  sailing,  there  are  deviations  from  the 
above  possibilities.  When  K  is  positive  and  the  sail  is  flipped  off,  FORCE  = 
CONST  and  CHAN  is  called.  When  a  simple  sail  is  used  FLAM  =  CONST  and 
return.  For  this  sail  there  is  a  constant  force  along  A.  When  a  fixed  sail  is 
used,  13  =  10,  FLAM  =  CONST  and  TACK  or  CHAN  is  called  depending  on  16. 
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Tack  calculates  the  forces  in  other  directions  as  in  addition  to  the  force 
along  X  the  FLAM  returns.  Control  returns  to  the  Main  program. 

When  K  is  negative  (or  zero)  the  force  FLAM  is  calculated  including  the 
CUTDOWN  or  not.  The  first  time  through  CON  =  CONST,  XLAM  is  corre¬ 
lated  by  CORR(TEL).  It  is  then  determined  if  the  modified  XLAM  lies  within 
the  band,  if  so  FLAM  =  0  and  control  is  returned  to  the  main  program.  If  the 
true  longitude  is  outside  the  desired  stationkeeping  range  control  passes  to 
statement  1  or  11  within  FLAM.  If  the  satellite  is  too  far  east,  FORCE  = 

-  CON,  the  dial  can  be  flipped  with  the  orbit  if  that  scheme  is  used,  or  the 
sail  can  be  offset  -90°  if  that  scheme  is  used.  If  the  satellite  is  too  far 
west,  FORCE  =  CON,  the  sail  can  be  flipped  against  the  orbit,  OFFSET  can 
be  0°.  In  either  case  control  passes  to  the  same  point  4,  where  K  will  re¬ 
turn  control  for  the  rest  of  the  orbit  as  long  as  the  impulses  are  being  given. 

At  statement  4  M  dutifully  counts  the  impulses.  CHAN  or  TACK  or 
neither  is  called  according  to  16.  If  M  exceeds  NUM,  the  fuel  consumption 
report  is  issued.  On  all  the  sail  schemes  the  fuel  is  set  to  zero.  And  the 
counters  are  reset  for  the  next  orbit  as  in  the  non-sailing  schemes. 

The  double  precision  function  FLAM(XLAM)  returns  to  the  main  program 
NUM  times  an  orbit  the  force  along  X  necessary  to  stationkeep  within  a  spe¬ 
cified  range  according  to  various  force  schemes  using  thrusters  and  solar 
sails . 
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APPENDIX  E 


SUBROUTINE  NEWFOR 


Subroutine  NEWFOR  is  an  alternative  method  of  calculating  the  force  nec¬ 
essary  to  keep  a  satellite  in  a  synchronous  orbit  using  the  bilateral  thrusting 
method  proposed  by  A.  Braga-Illa. 

A.  Theory 

Figure  18  shows  a  sample  of  longitude  vs.  time,  and  the  regions  ±  A,  ±B, 
±  C.  The  value  of  the  force,  con,  at  any  point  will  be  determined  by  which  re¬ 
gion  we  are  in  at  that  point  and  where  we  were  before. 

I  In  region  C  and  -C,  con  =  and  - C ^  resp. 

In  region  B  and  -B,  con  =  and  -C^>  resp. 

For  regions  A  and  -A  there  is  a  complicated  scheme  for  determining 
the  value: 

II  If  we  get  into  region  +A  via  +B  or  -  A  via  -B,  then  con  =  O. 

If  -f-A  via  -A,  con  =  or  -C^  resp. 

III  As  soon  as  we  get  the  sequence  (A-A+A)  or  (-A+A-A)  we  start 
damping  the  force  basically  as  follows.  Assume  we  get  the  sequence  +A-A+A 
and  then  the  curve  never  leaves  the  +A,  -A  region. 


Region: 

A 

B 

+A 

-A 

+A 

-A 

+A 

-A 

+A 

-A 

+A 

-A 

Force: 

C1 

C2 

0 

"C1 

C1 

2 

"C1 

2 

"C1 

4 

"C1 

4 

C1 

8 

"C1 

8 

C1 

TT~ 

-C1 

IT 

Region: 

+A 

-A 

+A 

-A 

+A 

-A 

Force: 

C1 

32 

-C1 

32 

0 

"C1 

32 

0 

"C1 

32 

where  A  is  the  region  where  we  start  damping.  After  (*)  we  do  not  damp  but 
keep  ±  C  ^/32,  0  as  long  as  we  stay  in  +A-A  and  do  not  stray  into  B. 
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TIME - ► 


Fig.  18.  Longitude  vs  time. 
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The  most  practical  way  to  show  the  deviations  from  the  basic  scheme  is 
via  figure  19-  Note  the  following: 

(1)  When  we  get  into  region  C  (see  arrow  1)  we  start  over  in  the  sense 
that  no  damping  occurs  until  we  get  the  ±A^A±A  sequence  again  (see  arrow  2). 

(2)  When  damping  has  started  and  we  go  into  region  ±B,  and  then  back 
into  ±A  we  keep  the  damped  value  that  we  had  before  going  into  region  ±B,  with 
the  appropriate  sign  (see  arrows  3,4). 

(3)  When  we  reach  the  maximum  damping  stage  ( ±0^/32,  0)  and  then 
go  into  region  ±B,  upon  re-entry  into  ±A  the  force  =  0  by  part  II  and  if  we  then 
go  into  T  A  from  there,  force  =  C  ^/32  with  the  appropriate  sign  ( see  arrows  5,  6) . 

B.  Data  Cards 

To  specify  that  subroutine  NEWFOR  is  to  be  used  to  calculate  the  force  we 
must  specify  four  quantities  on  the  data  cards: 

NEWFR  =  T,  Cj  =  ,  C2  =  ,  C3  = 

where  C^,  C^  are  the  values  assigned  to  the  force  in  regions  B  and  C  respect¬ 
ively  and  Cj  is  the  non- zero  undamped  value  of  the  force  in  region  A. 
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Fig.  19.  Example  of  possible  use  of  force  scheme  of  A.  Braga-Illa. 
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APPENDIX  F 


DOUBLE  PRECISION  FUNCTION  CORR  (TEL)  returns  to  FLAM  an  angle 
by  which  the  longitude  XLAM  is  thereby  modified.  In  FLAM:  X  =  CORR(TEL), 
XTEMP  =  XLAM  -  X.  The  program  calculates  the  apparent  right  ascension 
of  the  sun  and  the  time  of  ephemeris  transit,  good  for  the  year  1967,  and  accu¬ 
rate  to  about  four  seconds  of  arc  of  the  right  ascension  and  about  1.  5  seconds 
on  the  ephemeris  transit.  There  are  eight  options  in  this  program  determined 
by  the  13  variable.  Two  of  the  options  call  the  SUBROUTINE  CURT,  and  two 
of  the  options  call  the  SUBROUTINE  ALVISE.  When  13  =  4,  the  logic  in  FLAM 
bypasses  the  call  to  CORR. 

If  13  is  7  or  8,  immediately  after  converting  TEL  from  seconds  of  days, 
the  program  branches  to  110.  Day  is  modulated  by  365,  and  CURT  is  called 
to  determine  VALUE.  If  13  =  8,  VALUE  is  changed  to  radians  and  a  random 
error  is  added  on.  If  13  =  7,  CORR  is  simply  VALUE  in  radians  and  control 
is  returned  to  FLAM.  CURT  picks  up  an  ephemeris  sun  sensor  value  from 
the  step  function  tables  stored  in  the  BLOCK  DATA/SPETS/.  These  tables 
simulate  the  logic  in  an  electronic  sun  sensor  within  the  satellite. 

For  the  other  options,  the  program  calculates  the  mean  anomaly,  sun's 
coordinates,  referenced  to  a  geocentric  coordinate  system,  right  ascension 
of  the  sun,  right  ascension  of  the  mean  sun,  and  DIFF  the  difference  between 
the  right  ascension  of  the  sun  and  the  right  ascension  of  the  mean  sun  added 
to  6.  28  when  |DIFF  |  is  less  than  3.  Now  the  program  branches  according  to 
13  to  calculate  VALUE.  In  each  case  after  VALUE  is  calculated,  control  goes 
to  statement  31  where  CORR  =  DIFF  -  VALUE. 

If  13  is  5  or  6,  VALUE  is  returned  from  a  call  to  SUBROUTINE  ALVISE. 
When  13  =  6  DIFF  has  a  random  error  added  in.  ALVISE  as  a  subroutine  op¬ 
erates  identical  to  CURT,  looking  up  a  value  in  a  BLOCK  DATA.  When  con¬ 
trol  is  returned  to  CORR,  CORR  branches  to  31. 

If  13  =  1,  CORR  =  DIFF  and  control  is  returned  to  FLAM.  When  13  =  2, 
DAY  is  modulated  by  365.  The  formula  for  VALUE  is  determined  by  which 
part  of  the  year  DAY  falls. 
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Value 


Day 

0-45 
46  -  129 
130  -  211 
212  -  310 
310  -  360 


100  VALUE  =  .  08D0  *  (DAY  +  2.  2  5D0)  *  XX 
102  VALUE  =  -  .  0584D0  *  (DAY  -  111.  DO)  *XX 
104  VALUE  =  .  033D0*  (DA Y  -  1 62 .  DO)  *  XX 

106  VALUE  =  -0.  0584D0  *  (DAY  -  240.  DO)  *  XX 

107  VALUE  =  .  08D0  *  (DAY  -  363.  DO)  *  XX 


After  each  calculation  of  VALUE  control  branches  to  31.  When  13  is  3,  DIFF 
is  modified  using  the  random  error  function,  then  logic  follows  the  same  for 
13  =  2.  There  is  a  loop  for  13  =  4,  although  CORR  is  not  called  when  13  =  4, 
to  avoid  an  error  message  from  the  compiler. 
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APPENDIX  G 


The  SUBROUTINE  TACK  (XLAM)  and  its  ENTRY  CHAN  (XLAM)  return 
to  MAIN  via  FLAM  three  forces,  FORCCT,  FORLAM,  and  ALONGR  ,  which 
are  used  to  calculate  D2THET,  D2XLAM,  and  D2R  respectively.  The  ENTRY 
CHAN(XLAM)  fixes  the  sail  along  the  radius  to  the  satellite  from  the  earth. 
According  to  the  18  variables,  the  sail  remains  upright  on  the  radius  or  is 
flipped  throughout  part  of  the  orbit.  FORCCT  is  the  solar  force  exerted  by 
the  sail  in  the  direction  of  theta.  FORLAM  is  the  solar  force  exerted  by  the 
sail  in  the  X  direction.  ALONGR  is  the  solar  radiation  force  exerted  along 
the  radius.  L.  Black  has  pointed  out  that  the  most  efficient  sail  would  be  one 
of  high  thermal  emissivity  on  both  sides  with  a  high  solar  absorptivity  on  one 
side  and  a  high  solar  reflectivity  on  the  other.  A  sail  made  out  of  alumi¬ 
nized  microsheet  stuck  on  an  aluminum  sheet  painted  black  on  the  back  side 
would  be  a  good  approximation  to  this.  The  average  force  along  the  orbit  for 
a  sail  whose  plane  is  coincident  with  the  orbit  radius  vector  and  the  north- 
south  direction  is  given  approximately  by 


F 


AS 

2C 


where 

A  =  is  the  area  of  the  sail 
S  =  solar  constant 

C  =  velocity  of  light 


The  sail  placed  on  angle  bisector  of  angle  formed  by  sun  and  orbital  tangent 
in  direction  of  the  motion  is  the  most  efficient  sail;  the  program  allows  also 
an  offset  from  this  position,  OFFSET. 

When  the  solar  rays  hit  the  sail  they  generate  two  forces  which  must  be 
resolved  into  the  three  components  needed.  The  two  forces  are  one  a  force 
perpendicular  to  the  sail  itself  and  opposite  to  the  side  of  the  sail  hit  by  the 
solar  radiation;  the  second  is  a  force  along  the  sail  due  to  its  various 
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absorptivity  properties.  Both  sides  of  the  sail  are  reflective  when  the  sail  is 
"tacked"  and  when  it  is  flipped.  When  ENTRY  CHAN(XLAM)  is  used,  the  sail 
is  fixed  on  the  radius  from  the  earth  to  the  satellite,  and  one  side  is  highly  re¬ 
flective  while  the  other  side  is  highly  absorptive.  The  three  components  of  the 
generated  forces  are  the  force  along  the  orbit,  FORLAM;  the  force  along  the 
radius  vector  of  the  satellite,  ALONGR,  and  the  force  along  X,  FORLAM. 

In  this  subroutine  all  the  vectors  are  unit  vectors  and  lie  in  the  orbital 
plane  of  the  satellite.  The  coordinate  system  is  geocentric,  with  the  X-axis 
lying  along  Aries.  XNORMA,  YNORMA  are  the  coordinates  of  the  normal  to 
the  sail  itself.  The  coordinates  of  the  sun  are  calculated  as  in  Moulton  pp  171, 
and  then  the  vector  is  projected  into  the  orbital  plane  with  coordinates 
XPROJ,  YPROJ.  XFORCE,  YFORCE  are  the  coordinates  of  the  orbital  tan¬ 
gent.  XSAIL,  YSAIL  are  the  coordinates  of  the  normal  to  the  sail.  XRESAL, 
YRESAL  are  the  coordinates  of  the  sail  itself.  The  angle  between  the  sail  and 
the  radius  is  ANCHOR.  DOTPRO  is  the  cosine  of  the  angle  between  the  normal 
to  the  sail  and  the  projection  of  the  sun.  The  absolute  value  of  the  cosines  is 
used  to  simplify  the  logic  in  determining  quadrants.  ROCK  is  the  absolute  value 
of  the  cosine  of  the  angle  between  the  sail  and  the  radius.  TRAP  is  the  absolute 
value  of  the  cosine  of  the  angle  between  the  sun  and  the  normal  to  the  sail.  When 
DOTPRO  is  less  than  zero  the  sail  is  receiving  the  solar  force  on  its  absorptive 
side.  POT  is  the  force  exerted  perpendicular  to  the  sail  due  to  the  solar  radia¬ 
tion  pressure.  POT  =  (2  -  A)  *  DOTPRO  **  2.  ALOSAL  is  the  force  along  the 
sail  due  to  a  non-perfect  sail. 

ALOSAL  =  -A  *  TRAP  *  (XSUN  *  SRESAL  +  YSUN  *  YRESAL) 

There  is  a  logical  sequence  to  test  whether  the  sail  is  "ON"  with  the  flipped 
sail  logic.  If  it  is  on  one  continues,  if  not  the  program  branches  to  40,  where  a 
sequence  of  equations  set  up  the  small  forces  when  the  sail  is  turned  "OFF". 
Four  variables  are  used  to  set  up  the  signs  of  the  forces  according  to  the  quad¬ 
rant  the  sail,  sun,  and  orbital  tangent  are  in.  They  are  COD,  FISH,  HERRING, 
and  SALT.  They  are  initially  +1  and  are  made  -1  when  necessary.  The  solar 


56 


force  is  exerted  opposite  to  the  normal  to  the  sail  when  the  angle  between  the 
sun's  projection  and  the  sail  is  between  zero  and  180;  the  solar  force  is  ex¬ 
erted  along  the  normal  to  the  sail  when  the  angle  between  the  sun's  projection 
and  the  sail  is  between  180°and  360?  The  coordinates  of  the  solar  force  are 
XSORAF,  YSORAF.  HERRING  is  negative  when  the  angle  between  the  solar 
force  and  the  orbital  tangent  is  greater  than  180.  SALT  is  negative  when  the 
angle  between  the  sail  and  the  orbital  tangent  is  greater  than  180.  COD  is 
negative  when  the  angle  between  the  solar  force  and  the  radius  is  greater  than 
180.  FISH  is  negative  when  the  angle  between  the  radius  and  the  sail  is  greater 
than  180.  The  three  components  of  the  solar  radiation  force  acting  on  the  sat¬ 
ellite  are: 

FORLAM  =  HERRING  *  POT  *  ROCK  -f  SALT  *  ALOSAL  *  BOAT  |  0 

I  bee 

ALONGR  =  COD  *  POT  *  BOAT  -f  FISH  *  ALOSAL  *  ROCK  Fig- 

FORCCT  =  - AB1  *  TRAP  *  ZSUN  "q6 

When  the  sail  is  flipped  to  its  "OFF”  position  there  are  slight  forces  acting 
upon  the  satellite.  DOTPR1  is  the  cosine  of  the  angle  between  the  sun  (not  its 
projection)  and  the  sail.  DOTPR2  is  the  absolute  value  of  ZSUN.  FORCCT 
is  negative  when  ZSUN  is  positive.  Otherwise  the  three  forces  are: 

FORLAM  =  AB1  *  DOTPR2  *  DOTPRO 
ALONGR  =  -  AB1  *  DOTPR2  *  DOTPR1 
FORCCT  =  (2  -  AB1)  *  DOTPR2 

AB1  is  the  absorptivity  received  by  the  program  as  A1  and  A2  and  set  equal 
to  the  proper  one  depending  on  the  quadrant  one  is  in.  Following  is  a  diagram 
of  the  forces  for  an  offset  sail  in  one  quadrant. 
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FORCE  EXERTED  BY  SUN  ALONG  NORMAL  TO  SAIL 


Fig.  20.  VECTOR  DIAGRAM  OF  FORCES  acting  on  solar  sail. 
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subroutine  talic  flow  cmart 
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APPENDIX  H 


SUBROUTINE  MINMAX  (X,  N,  XMIN,  XMAX) 


Subroutine  minmax  is  designed  to  find  the  maximum  and  minimum 
value  of  a  set  of  N  numbers. 


60 


61 


APPENDIX  I 


SUBROUTINES  ALVISE  AND  CURT 


SUBROUTINE  ALVISE  contains  the  information  for  the  step  function  ap¬ 
proximation  to  the  equation  of  time  used  for  the  LES-6  satellite. 

SUBROUTINE  CURT  contains  the  information  for  the  step  function  ap¬ 
proximation  to  the  equation  of  time  used  for  the  LES-6  satellite  analemma 
sensor . 
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5-  C.3-  &  AB9 
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3-43-6*10 
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appendix  j 


LISTING  OF  THE  PROGRAM 
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